Origin of the structural phase transition in Li 7 La 3 Zr 2 Oi2 
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Garnet-type Li7La3Zr20i2 (LLZO) is a solid electrolyte material with a low-conductivity tetrag- 
onal and a high-conductivity cubic phase. Using density-functional theory and variable cell shape 
molecular dynamics simulations, we show that the tetragonal phase stability is dependent on a 
simultaneous ordering of the Li ions on the Li sublattice and a volume-preserving tetragonal dis- 
tortion that relieves internal structural strain. Supervalent doping introduces vacancies into the Li 
sublattice, increasing the overall entropy and reducing the free energy gain from ordering, eventually 
stabilizing the cubic phase. We show that the critical temperature for cubic phase stability is low- 
ered as Li vacancy concentration (dopant level) is raised and that an activated hop of Li ions from 
one crystallographic site to another always accompanies the transition. By identifying the relevant 
mechanism and critical concentrations for achieving the high conductivity phase, this work shows 
how targeted synthesis could be used to improve electrolytic performance. 



The garnet structure material Li7La3Zr20i2 (LLZO) 
has an ionic conductivity that varies by two orders 
of magnitude depending on whether synthesis pro- 
duces a Cubic ((7cubic = 1.9xl0" 4 S/cm) [H' or tetragonal 

(at ctra — 

1.63xl(T 6 S/cm) Q phase. With the higher con- 
ductivity, LLZO has excellent solid electrolyte material 
characteristics. Its stability against both Li metal and 
standard cathode LiCo02 |3|, combined with a ~5 eV 
band gap [l[ and high ionic conductivity, make it suitable 
for exploiting the full voltage difference between anode 
and cathode while circumventing the safety concerns in- 
herent to all current liquid electrolytes. Initially, the pro- 
cess by which the cubic ground state could be stabilized 
against the tetragonal was unknown and uncontrolled. 
The relevant factor was eventually traced down to the 
addition of Al in the structure, whether via accidental 
uptake from Al-containing crucibles Q or via intentional 
doping 0, Inclusion of other supervalent ions, in- 
cluding Ta, Nb, and Ga @, 0], has also been successful 
in producing the high conductivity phase. However, the 
underlying mechanism that controls the transition has 
so far remained a mystery, hindering further progress to- 
ward improving the material for practical usage. 

Here we use our recently developed variable cell shape 
density-functional theory (DFT) plus molecular dynam- 
ics (MD) method to investigate the driving force behind 
the tetragonal to cubic transition which subsequently 
raises the conductivity. We find that in the cubic phase, 
the Li sublattice is always disordered (all Li symmetry 
sites partially occupied), while in the tetragonal phase it 
is always ordered (all Li sites either full or empty). We 
find that the energy is lowered by a simultaneous order- 
ing of the Li atoms that relieves Li— Li Coulomb repulsion 
but unfavorably distorts the ZrOg octahedra and a lattice 
distortion that restores the preferred Zr— O bonds. The 
two symmetry-breaking but volume-preserving phenom- 



ena always occur in conjunction and either one alone ac- 
tually raises the total energy of the system. When Al 3+ is 
doped into the system, charge compensation takes place 
through creation of Li + vacancies that reduce the free en- 
ergy advantage of complete ordering on the Li sublattice, 
eventually leading to disorder and a transition to cubic 
symmetry. The transition is always signaled by a sudden 
shift of Li occupation which can be used to pinpoint the 
critical temperature and vacancy concentration. Using 
these criteria, we estimate that the critical concentration 
of Al dopants necessary to achieve the high-conductivity 
cubic phase of Li7_2 a: AL E La3Zr20i2 is x = 0.2 (vacancy 
number n vac = 2x — 0.4 per formula unit) at zero tem- 
perature, in good agreement with experiment. We show 
that the cubic phase can be reached at some tempera- 
ture, regardless of Li content, but that the critical tem- 
perature drops as a function of vacancy number. The 
understanding uncovered in this work will be useful for 
choosing more efficient dopants and further raising the 
ionic conductivity. 

The Li sublattices in the cubic and tetragonal LLZO 
phases are shown in Fig. [TJ Li positions in each struc- 
ture are generally referred to as Li(l) if they are tetra- 
hedrally coordinated to oxygen, and Li(2) and Li(3) if 
they are octahedrally coordinated. To avoid confusion, 
we refer to each site using its crystallographic notation, 
with superscript c or t to designate cubic and tetrag- 
onal lattices, respectively. Each tetrahedral cubic 24d c 
site is surrounded by four pairs of octahedrally coordi- 
nated 96h c sites, and all sites are partially occupied (see 
Fig. [1]). Because of Coulomb repulsion, it is energetically 
prohibited for both members of each pair of 96/i c sites 
to be occupied, and if a particular 24<i c site is occupied, 
the adjacent 96/i c sites are consistently unoccupied [||. 
The tetragonal distortion transforms the cubic 24d c sites 
into fully occupied 8a* sites, often denoted as Li(l), and 
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FIG. 1: (Color online) Li sublattice in the cubic (left) and 
tetragonal (right) phases of LLZO. All Li positions are in- 
cluded, although in the cubic phase not all are occupied. The 
Li(l) atoms (8a* and 24d c ) are large gray (gold online), Li(2) 
or 16/* are white, and Li(3) or 32g % are dark gray. The cubic 
Li(l) positions that become vacant upon transition to the or- 
dered tetragonal structure (16e t ) are indicated by small (gold 
online) spheres. 



unoccupied 16e* sites. The 96h c sites are transformed 
into fully occupied 16/' and 32g t sites, often denoted as 
Li(2) and Li(3), respectively, with the rest fully unoccu- 
pied. As we will show, a shift from 8a* sites to 16e* sites, 
both subsets of 24d c , always accompanies the tetragonal 
to cubic transition. 

All of our simulations use density-functional theory, 
with the Perdew-Burke-Ernzerhoff exchange-correlation 
functional 0, [lfl ] as implemented in the Vienna Ab-initio 
Simulation Package (VASP) projector-augmented- wave 
software version 5.2.12 [ll],[l2j]. The calculations used an 
energy cutoff of 400 eV, except for the relaxations used to 
explain the microscopic basis of the structural transition, 
which used a cutoff of 600 eV. Molecular dynamics (MD) 
time evolution is carried out using the velocity Verlet al- 
gorithm, as implemented in the libAtoms [131 ] package, 
with a time step of 1 fs and using forces and stresses cal- 
culated at each time step with VASP. All MD simulations 
use constant temperature and stress. The time propaga- 
tion algorithm, specified in Supplementary Information, 
is a rediscretization, using the ideas from Ref. [LJ and 



Ref. Il5. of the Langevin constant pressure algorithm in 




for a modified version of the equations of motion 
171 Each simulation starts from a unit cell of the 
experimental tetragonal structure of Ref. 0, for a total of 
8 formula units. The structure is relaxed with respect to 
ionic positions and unit cell size and shape using VASP, 
and then simulated at finite temperature and zero stress. 
While the MD simulation results implicitly include the 
effects of entropy, all energies we calculate explicitly and 
quote here include only DFT total internal energy. The 
system is allowed to evolve for at least 24 ps with ionic 
motion as well as cell shape and volume changes, so it 
can spontaneously transform into other structures. To 
determine the effect of composition and temperature we 



simulate the stoichiometric system, with 56 Li atoms, as 
well as systems with 1, 2, and 4 Li + ion vacancies per 
simulated cell, at temperatures ranging from 300 K to 
1300 K. The net negative charge of the vacancies is com- 
pensated by a uniform positive background charge. The 
vacancies are formed by removing Li + ions randomly se- 
lected from the tetragonal 16/' and 32g t sites, where the 
vacancy formation energy is about 100 meV lower than 
at the 8a* sites. The crystallographic site identity of each 
atom is determined during each MD trajectory (see Sup- 
plementary Information). 

An example of the time evolution of one system, with 
n vac = 0.25 and T = 600 K, is plotted in Fig. H The ra- 
tios of the lattice constants along x and y to that along z 
initially show the tetragonal structure with one axis (a z ) 
about 3% shorter than the others. The system trans- 
forms into a cubic phase where both axis ratios fluctuate 
around 1 at t ~ 5 ps. It fluctuates back to a tetragonal 
phase at t ~ 15 ps, and then again to cubic at t ~ 30 ps 
where it remains for the rest of the simulation. The vol- 
ume is not affected by the unit cell shape change. We 
also plot the occupancies of various crystallographic sites 
of the cubic and tetragonal structures, computed as de- 
scribed above. We see a perfect correlation between the 
symmetry of the unit cell parameters and the occupancy 
of the tetragonal sites. The 8a* occupancy is initially 
near 1, as it is in the experimental structure, and most 
of the remaining atoms are identified as 16/' and 32g t . 
When projected into the cubic structure the 8a* atoms 
are identified as 24<i c and the 16/* and 32g t atoms are 
identified as 96ft°, as expected by symmetry. Whenever 
the system transforms into the tetragonal structure the 
occupancies of 8a' and 16/ t +32g t sites drop. Since we do 
not see a corresponding change in 24d c and 96/i c , we con- 
clude that the atoms are moving from the subset of the 
cubic sites that are occupied in the tetragonal ordering 
to the other cubic structure sites. This is indeed seen in 
the increased occupancy of 16e' sites, which correspond 
to 2Ad c sites that are unoccupied in the experimental 
tetragonal structure. 

The unit cell shape, 8a* occupancy, and unit cell vol- 
ume, averaged over the last 5 ps of each run, are listed 
m Table H At each vacancy concentration the system 
transforms from an ordered tetragonal structure at low 
T to a cubic structure at higher T, and the transition 
temperature goes down with increasing vacancy concen- 
tration. Note that at n vac = 0.5 and T — 300 K our sim- 
ulations predict an ordered tetragonal structure that is 
different from the experimentally observed low T tetrag- 
onal structure for the stoichiometric composition. The 
structure we find has one long and two short axes, an 
occupation of the original tetragonal 8a* sites of 0.5, and 
an occupation of 1 of the original tetragonal 16e' sites 
(which are unoccupied in the stoichiometric tetragonal 
structure). The volume is much more sensitive to tem- 
perature than to vacancy number, with a linear thermal 



3 




5 10 15 20 25 30 35 40 45 50 
t(ps) 

FIG. 2: (Color online) Evolution over time of structural and 
site occupation quantities for a sample system with n V ac = 2 
at T — 600 K. Top: unit cell shape (a x /a z blue, a y /a z red) 
and volume (black). Middle: 96h c (black) and 16/ t +32p t 
(red) lattice site occupation. Bottom: 24d c (black), 8a* (red) 
and 16e* (blue with symbols) lattice site occupations. 

expansion coefficient (from a linear fit of V 1 ^ 3 vs. T to 
the n V ac = 0.5 results, where we have the widest range 
of temperatures) of 2.2 x 1CT 4 K" 1 . The vacancy num- 
ber dependence (from a linear fit of V vs. n vac to the 
T = 600 K and T = 800 K results), which is the va- 
cancy formation volume, is 10±1.5 A 3 /vacancy. Note 
that this value is dependent on the charge compensation 
mechanism, which is a uniform background charge in this 
calculation, compared to any of a variety of supervalent 
dopants in experiment. In all cases, a transition to cu- 
bic symmetry is preceded by a sudden decrease in the 
8a* occupation and any return to tetragonal symmetry 
is characterized by a reoccupation of these sites. 

Wc have also computed the energy difference between 
the tetragonal and cubic structures as a function of va- 
cancy concentration. For the tetragonal structure we re- 
laxed 10 structures with random vacancies at 16/* and 
32g* sites. For the cubic structure, which is disordered 
even at the stoichiometric composition, we use 10 con- 
figurations from a uniform sampling of all configurations 
that obey the restriction that no pairs of nearest neigh- 
bors (adjacent 96h° pairs, or a 24d c and its nearest neigh- 
bor 96h c ) are simultaneously occupied. We plot the en- 
ergy difference between the minima and means of each 
of the ten structures in Fig. [3l The tetragonal structure 
is energetically favored for n vac < 0.25, and the cubic 
is favored for n vac > 0.5. Entropy effects will shift the 
equivalent curves for the free energies. We expect that 
this shift will reduce the free energy of the tetragonal 
phase more than the cubic phase, because the tetrago- 
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FIG. 3: (Color online) Energy difference between tetragonal 
and cubic structures as a function of vacancy number, for 
minimum energy configuration (solid red line) and mean con- 
figuration energy (dashed blue line). 
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FIG. 4: Relaxed energies (closed symbols) of 10 disordered 
configurations and one at the experimentally observed tetrag- 
onal phase order, for a cubic unit cell and for a fully relaxed 
one, relative to the lowest energy fully relaxed ordered tetrag- 
onal cell. The energy of a model including only ZrOe distor- 
tions is also plotted (open symbol). 

nal phase starts out ordered and therefore gains more 
entropy with the introduction of vacancies than the al- 
ready disordered cubic phase. The transition vacancy 
concentration will therefore shift to higher values as the 
temperature increases, although we expect this shift to 
be small at T = 300 K. 

To understand the relationship between Li or- 
der/disorder and the tetragonal/cubic lattice parameters, 
we perform total energy calculations for the stoichiomet- 
ric system in 10 disordered configurations (again from a 
uniform sampling of configurations obeying first-neighbor 
exclusion) and in the experimentally observed Li order in 
the tetragonal phase. We first relax each one with respect 
to atomic positions and unit cell volume but constrain the 
unit cell shape to remain cubic, and then relax with full 
freedom for the unit cell size and shape as well as atomic 
positions. The total energies are plotted in Fig. @] The 
energies of the disordered system are generally highest, 
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TABLE I: Mean unit cell distortion (a a /a z — 1 for a—x and y, in %) (top) and occupancy of 8a* sites per formula unit (bottom), 
averaged over the last 5 ps of each run, for different vacancy numbers n va c and temperatures. Transition temperature T c is also 
indicated for each quantity with a clear signal of a structural transition (see text). 
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with only a small gain (and correspondingly small change 
in unit cell shape) from relaxing the cubic unit cell con- 
straint. Imposing the Li order while maintaining a cubic 
unit cell actually increases the energy slightly as com- 
pared with the mean disordered energy. Only once the 
ordered cell is allowed to relax to the tetragonal shape 
does its energy become lower than the disordered system, 
and the ordered system emerges as the T — K ground 
state. 

To understand the source of the coupling between en- 
ergy and structure, we calculated the pair distribution 
functions (PDF) of various ion pairs in the system for 
the cubic and tetragonal ordered systems, along with a 
simple point charge model, based on the nominal ionic 
charges. The point charge model shows that the overall 
Coulombic energy is actually increased upon relaxation 
to a tetragonal unit cell, indicating that the energy gain 
is elsewhere. The PDFs show that the La 3+ — 2 ~ dis- 
tances do not vary, presumably since these are the two 
largest ionic charges in the system. Li— O and Li— Li dis- 
tances do change, but the highly ionic nature of these 
interactions means that such changes are nearly entirely 
accounted for in the point charge model. The shifting of 
Li— Li and Li— O distances, combined with the rigidity 
of the La— O spacing, results in some distortion of the 
ZrOg octahedra. Unlike the other pairs, the Zr— O in- 
teraction is at least partially covalent. In the cubic cell, 
the Zr-0 bond lengths are 2.130±0.020 A and O-Zr-0 
bond angles are 180±4.0°. The relaxation to a tetrag- 
onal cell relieves this distortion, restoring the octahe- 
dra to a uniform Zr-0 bond length of 2.125±0.005 A 
and a O-Zr-0 bond angle of 180±0.01°. To estimate 
the energetic contribution of ZrOe octahedra distortions 
we parametrize the calculated total energies of the or- 
dered tetragonal cell as a quadratic function of Zr— O 
bond length and O— Zr— O bond angle by computing the 
energies for small displacements of an O atom. This 



parametrization predicts an energy for the ordered cubic 
structure of 0.083 eV/formula unit relative to the tetrag- 
onal structure (Fig. @|, very nearly equal to the DFT 
energy difference. We therefore conclude that lithium 
ordering to relieve Li-Li Coulomb repulsion leads to in- 
ternal rearrangements of ions to maintain La 3+ — 2+ dis- 
tances, which leads to a lattice distortion that allows the 
ZrOg octahedra to preserve their preferred shape. 

In summary, our variable cell shape DFT MD sim- 
ulations of LLZO show that the DFT ground state at 
low temperatures is the experimental ordered tetragonal 
structure, and that at higher temperatures the system 
transforms into the experimental disordered cubic struc- 
ture. The transition temperature decreases with increas- 
ing Li + vacancy concentration, and the disordered cubic 
structure has lower energy when the number of vacancies 
per formula unit is larger than about 0.4. These results 
are in agreement with the experimental observations of a 
transition as a function of Li + vacancy concentration Q . 
The microscopic cause of the structural transition is the 
coupling between the unit cell shape and the hopping of 
atoms from the subset of the disordered cubic sites occu- 
pied in the ordered tetragonal structure to the remain- 
ing cubic sites, which are unoccupied in the tetragonal 
structure. The relative stability of the ordered tetrago- 
nal low temperature structure is driven by the ordering, 
which reduces Li— Li Coulomb repulsion but distorts the 
ZrOe octahedra, and the tetragonal distortion which al- 
lows these octahedra to return to their preferred high- 
symmetry shape. Our simulations enable us to explain 
the atomistic mechanism behind this finite temperature 
structural transformation in a complex material, and pre- 
dict the number of vacancies necessary to achieve the 
high conductivity material. This should enable better 
doping schemes that optimize ionic conductivity by pro- 
viding the requisite number of vacancies with as few arti- 
ficial dopants as possible, thereby realizing the potential 
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of LLZO as a truly stable solid electrolyte material. 

This work was supported by the Naval Research Labo- 
ratory core 6.1 research program and the Nanoscience In- 
stitute. Computer time was provided through the DOD 
HPCMPO at the ERDC and AFRL MSRCs. 



* Resident at the Center for Computational Materials 
Science, Naval Research Laboratory, Washington, DC 
20375. 

[1] R. Murugan, V. Thangadurai, and W. Weppner, Angew 

Chem Int Edit 46, 7778 (2007). 
[2] J. Awaka, N. Kijima, H. Hayakawa, and J. Akimoto, J. 

Solid State Chem. 182, 2046 (2009). 
[3] S. Ohta, T. Kobayashi, and T. Asaoka, J Power Sources 

196, 3342 (2011). 
[4] C. A. Geiger, E. Alekseev, B. Lazic, M. Fisch, T. Arm- 

bruster, R. Langner, M. Fechtelkord, N. Kim, T. Pettke, 

and W. Weppner, Inorg Chem 50, 1089 (2011). 
[5] Y. Jin and P. J. McGinn, J Power Sources 196, 8683 

(2011). 



[6] J. L. Allen, J. Wolfenstine, E. Rangasamy, and 
J. Sakamoto, J Power Sources 206, 315 (2012). 

[7] S. Adams and R. P. Rao, J Mater Chem 22, 1426 (2011). 

[8] H. Xie, J. A. Alonso, Y. Li, M. T. Fernandez-Diaz, , and 
J. Goodenough, Chem. Mat. 23, 3587 (2011). 

[9] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys Rev Lett 

77, 3865 (1996). 

[10] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys Rev Lett 

78, 1396 (1997). 

[11] G. Kresse and J. Hafner, Phys Rev B 47, 558 (1993). 
[12] G. Kresse and J. Furthmuller, Phys Rev B 54, 11169 
(1996). 

[13] http://libatoms.org/. 

[14] A. Jones and B. Leimkuhler, J Chem Phys 135, 084125 
(2011). 

[15] B. Leimkuhler, E. Noorizadeh, and O. Penrose, J Stat 

Phys 143, 921 (2011). 
[16] D. Quigley and M. I. J. Probert, J Chem Phys 120, 11432 

(2004). 

[17] E. B. Tadmor and R. E. Miller, in Modeling Materi- 
als (Cambridge University Press, Cambridge, 2011), pp. 
520-533. 



